Zfp335 establishes eTreg lineage and neonatal immune tolerance by targeting Hadha-mediated fatty acid oxidation

Regulatory T cells (Tregs) are instrumental in maintaining immune tolerance and preventing destructive autoimmunity, but how heterogeneous Treg populations are established remains largely unknown. Here, we show that Zfp335 deletion in Tregs failed to differentiate into effector Tregs (eTregs) and lose Treg-suppressive function and that KO mice exhibited early-onset lethal autoimmune inflammation with unrestricted activation of conventional T cells. Single-cell RNA-Seq analyses revealed that Zfp335-deficient Tregs lacked a eTreg population and showed dramatic accumulation of a dysfunctional Treg subset. Mechanistically, Zfp335-deficient Tregs displayed reduced oxidative phosphorylation and dysfunctional mitochondrial activity. Further studies revealed that Zfp335 controlled eTreg differentiation by regulating fatty acid oxidation (FAO) through direct targeting of the FAO enzyme Hadha. Importantly, we demonstrate a positive correlation between ZNF335 and HADHA expression in human eTregs. Our findings reveal that Zfp335 controls FAO-driven eTreg differentiation to establish immune tolerance.


Flow cytometry
For analysis of surface markers, cells were stained in PBS containing 2% FBS on ice for 30min.For analysis of intracellular cytokine staining, cells were stimulated for 4 h in vitro with PMA/Ionomycin in the presence of brefeldin A and monensin.The stimulated cells were fixed and permeabilized using a Fixation/Permeabilization Solution Kit (Biolegend).
Treg cells were incubated with 3uM BODIPY-FL-C16 and 2 μM BODIPY-493/503 for 30 min at 37℃.After the incubation, the cells were washed in RPMI 1640 medium with 2% FBS and continued for staining surface markers at room temperature in dark.
All samples were analyzed using a CytoFLEX flow cytometer (BECKMAN COULTER).FlowJo software and CytExpert software were used for data analysis.Cells were sorted with a FACSAria (BD Biosciences).

RNA isolation and quantitative PCR (qPCR)
Total RNA was extracted with a RNeasy Mini Kit (Qiagen), according to the manufacturer's instruction.cDNA was reverse transcribed using the cDNA synthesis kit (TOYOBO) and amplified with SYBR Green RT-qPCR Mastermix (GenStar) at StepOnePlusTM Real-Time PCR System (ThermoFisher).Primer sequences used in this study were summarized in Supplemental table 5.

HE staining and histopathology
Lungs, livers, kidneys, ears, pancreas, salivary glands, thymi, hearts and lacrimal glands were removed from 3-week-old WT and KO mice.Colons were removed from 3-week-old WT and KO mice or mice as mentioned in colitis model.Samples were formalin fixed, paraffin embedded and stained with haematoxylin and eosin before tissue histology.Photomicrographs were taken at x20 or x5 magnifications.

ELISA for autoantibody
Serum samples were collected from 3-week-old WT and KO mice.
Autoantibodies (anti-dsDNA) were measured using a detection kit from Alpha Diagnostic International (5110) according to the manufacturer's instructions.

Malate supplementation in vitro
For malate treatment in vitro, resting Treg (rTreg) cells purified from Tamoxifen treated-ER Cre and ER Cre Zfp335 fl/fl mice were sorted and activated with 5μg/ml anti-CD3 Ab, 2μg/ml anti-CD28 Ab and 500U/ml IL-2, and supplemented with or without 30 mM malate.72h later, cells were collected and analyzed by flow cytometry.

In vitro fatty acid supplementation
Sodium oleate (Sigma) was dissolved in PBS and stocked at 25 mM.Oleate was then dissolved by heating in a metal bath at 70℃ and conjugated with RPMI 1640 medium supplemented with 1.6% FA-free BSA.Purified rTreg cells were activated by 5μg/ml anti-CD3 Ab, 2μg/ml anti-CD28 Ab and 500U/ml IL-2 for 2-3 days in the presence or absence of 50uM sodium oleate.The cells were then collected and applied for flow cytometry analysis.

Bulk RNA-seq data analysis
Total RNA was isolated form CD4 + CD25 + YFP + Treg cells of 2-week-old Foxp3 Cre and Foxp3 Cre Zfp335 fl/fl mice and used for RNA sequencing analysis.Firstly, the bulk RNA-seq data were filtered by using SAOPnuke (version 1.5.6)(1) with parameters "-l 15 -q 0.2 -n 0.05 -Q 2".After removing low-quality bases RNAs, the clean data were mapped to mouse genome (mm10) by using HISAT2 (2) with parameters "-k 1 -p 4 -q --no-unal --dta --un-conc-gz".Then the expression levels of each gene were calculated by the transcripts per kilobase of exon model per million mapped reads (TPM) by using StringTie (3) with parameters "-t -C -e -B -A -p 1".The final TPM matrix of all samples was used for subsequent analysis.A 1.5-fold variance in expression levels, a P value less than 0.05, and an adjusted P value less than 0.1 were used as cutoffs to define differentially expressed genes.The P value and adjusted P value were calculated using R software (DESeq2) (4).

Single-cell RNA sequencing processing
The spleens were dissociated into single-cell suspensions with the following procedure: Spleens were processed with the flat end of a syringe in a 100 mm culture dish containing 5 ml cold FACS buffer (2% FBS in PBS), then passed through a 70 μm cell strainer into a 15 ml tube.Cells were centrifuged to remove the supernatant.Cell pellets were treated with 1ml ACK (Ammonium-Chloride-Potassium) Lysing Buffer to remove the red blood cells.After washing with 10 ml cold FACS buffer, the remaining cells were stained with 7AAD (Part 76332; Lot B226294 Biolegend) for 30 min at 4 °C before flow cytometric sorting using FACS Aria II Cell Sorter (BD Biosciences).The sorted CD4 + YFP + 7AAD -cells with a viability higher than 90% were used for 10X genomics scRNA-seq.
Furthermore, the single-cell library preparation was constructed using 10X Chromium Single Cell V3 Reagent Kits according to the manufacturer's protocol.
Cell Ranger (V5.0.1, https://support.10xgenomics.com/)was used to process scRNA-seq data and generate the matrix data containing gene counts for each cell per sample.Briefly, the 10X sequencing data were mapped to the mouse genome (mm10) which downloaded from 10X Genomics and generated the unique molecular identifiers (UMI) matrix of each cell by using Cell Ranger (version 5.0.1)count pipeline.

GO, KEGG and Gene Set Enrichment Analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using clusterProfiler (V3.18) package using genes specifically expressed in indicated Treg cell cluster.Gene Set Enrichment Analysis (GSEA) analysis was performed for each cell subpopulation using the scaled gene expression matrix and GSEA package (V.

Single-cell trajectory analysis
To reveal the differentiation relationship of various Treg subsets, Monocle (v3) (5) was used for pseudotime analysis.The Seurat object was converted to a Monocle3 object using as.cell_datA_set function.Then cluster_Cells and learn_graph functions were used to construct developmental trajectories in UMAP.The get_earliest_principal_node helper function as was used to assign a node for which the highest fraction of closest cells belonged to the rTreg cluster as the root node.Then, the 'order_cells' function was used to order cells and the plot_cells function was used to visualize the trajectory in two-dimensional spaces.

Hallmark Gene Set score quantification
To score individual cells for Hallmark pathway activities, we used multi-previous described methods analyzing different Treg subsets data.Firstly, the mouse Hallmark Gene Sets transformed from the human genes were used from msigdbr package, and gene sets were then used to score each cell.To eliminate the bias of sample background information, we selected gene set enrichment analysis methods based on single cell gene expression ranking AUCell (6), UCell (7), singscore (8) and ssGSEA (9).Of note, ssGSEA cancels the final standardization step, making it closer to the gene set enrichment analysis of a single cell.In addition, to evaluate whether the gene set is enriched in a certain cell subpopulation, we calculated the differential gene set in the enrichment score matrix by Wilcox test (the filter criterion for differential genes is that the P value after correction is less than 0.05).Finally, we used the rank aggregation algorithm (RRA) in the RobustRankAggreg package (10) (version 1.1.0)to comprehensively evaluate the results of the difference analysis, and screen out the genes that are significantly enriched in most gene set enrichment analysis methods Set (the filter criterion for comprehensive evaluation is P value less than 0.05).

2.5×10 7
Treg cells were sorted from the lymph nodes of WT mice on FACSAria Ⅱ (BD Bioscience) cell sorter.Millipore 17-10085 Chromatin Immunoprecipitation (ChIP) kit and anti-Zfp335 antibody (Novus) were used in the ChIP assay.Immunoprecipitated DNA was used for Illumina ChIP-seq sample preparation.In brief, 2.5×10 7 cells were crosslinked to chromatin with 1% formaldehyde.Reaction was stopped with 0.125M glycine.The cells were resuspended in cold nuclear lysis buffer and the chromatin was sonicated to yield fragments of ~300-500bp size, followed by overnight incubation with immunoprecipitation-grade anti-Zfp335 antibody and Magnetic Protein A/G Beads.The following day, beads were sequentially washed by low-salt, highsalt, LiCl, and TE buffers.Bound complexes were eluted in 150μl of elution buffer at 62°C for 2h with shaking, followed by reversal of formaldehyde crosslinking at 95°C for 10 minutes.DNA was eventually purified with spin columns.Immunoprecipitated DNA concentration was detected by the Qubit DNA broad range assay in the Qubit Fluorometer (Invitrogen).10ng immunoprecipitated DNA was prepared for sequencing using the Illumina ChIP-seq sample preparation protocol.The library products were enriched quantified and finally sequenced on Novaseq 6000 sequencer (Illumina) with PE150 model.Raw sequencing data were filtered by Trimmomatic (version 0.36).FastQ reads were aligned to the ensemble mouse genome (GRCm38) with STAR software (version 2.5.3a) using default settings.The MACS2 software (Version 2.1.1)was used to process peak calling.The library products corresponding to 200-sequencer (Illumina) with PE150 model.Genomic graphs were generated and viewed with the IGV (Integrative Genomics Viewer).